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Abstract 


The  purpose  of  this  work  was  to  mathematically  model  the  heat  transfer 
in  an  active  packed  bed  heat  exchanger  and  to  develope  an  appropriate  num¬ 
erical  scheme  for  the  solution  of  the  temperature  field.  As  active  packed 
beds  are  being  considered  for  use  in  nuclear  reactors,  concern  is  generated 
as  to  whether  regions  of  relatively  high  temperatures  will  develope.  Cuch 
concern  is  the  impetus  behind  the  application  of  our  program  to  two  appli¬ 
cations  concerning  heat  transfer  in  packed  beds  assuming  (i)  variable  por¬ 
osity  of  the  Brosilow  type  and  assuming  (ii)  the  existence  of  blockages. 

The  results  show  that  for  the  case  of  a  packed  bed  with  the 
porosity  distribution  the  so-called  hot  spot  occurrs  near  the  wall  of  the 
bed.  Its  magnitude  and  position  were  seen  to  be  strongly  affected  by  the 
assumed  velocity  profile  and  wall  boundary  conditions.  In  addition  as  the 
proposed  mathematical  equations  include  axial  conduction  as  opposed  to  a 
previously  developed  program  which  excluded  axial  conduction  and  as  the 
results  of  the  two  programs  approached  one  another  as  the  Re  increased,  it 
is  reasonable  that  the  previous  program,  (as  it  is  based  on  parabolic,  not 
elliptic  equations  as  in  the  present  case),  could  be  implimentod  with 
little  error  expected  due  to  the  lack  of  axial  conduction. 

In  the  case  of  the  blockage,  however,  because  of  small  mass  flow 
in  the  blocked  region,  the  axial  conduction  is  important  and  the  present 
program  must  be  used.  Effects  of  variable  or  constant  eddy  diffusivity, 
wall  boundary  conditions,  position  of  blockage  and  effective  thermal  con¬ 
ductivity  of  the  solid  are  investigated. 
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.  Tntrodv.ct  5  on 

'eat  cxc  angers  'ave  diversified  vises  in  industry. 

’.cat  to  a  flowin'’  meduiitm;  the  storage  of  heat;  and  c  emical  regeneration 
are  a  few  such  uses.  More  specifically,  packed  bed  heat  exc  angers  offer 
the  advantage  that  large  rates  of  heat/chemical  transfer  can  take  place  in 
a  relatively  small  volume.  A  packed  bed  is  typically  a  collection  of  small 
spheres,  (however  various  shapes  are  possible),  housed  in  a  cylindrical 
tube  in  such  a  manner  that  the  spheres1  positions  are  stationary;  see 
figure  1.  Small,  tightly  packed  particles  are  what  yield  t he  high  surface 
to  volume  ratio,  6*(l-fc)/Dp  ;  thereby  resulting  in  good  heat  transfer 
properties.  Such  an  advantage  is  why  packed  beds,  with  elements  coated 
with  fissionable  material;  Were  considered  as  likely  candidates  for  use 
as  nuclear  reactor  cores.  In  fact,  plants  using  such  packed  beds  as  reactor 
cores,  are  in  operation  presently. 

One  concern,  however,  with  the  application  of  packed  beds  in  nuclear 
reactors,  is  the  existence  of  regions  where  temperatures  can  achieve  sig¬ 
nificantly  higher  levels  than  tie  rest  of  the  field.  Such  concern  stems 
from  the  fear  that  greater-than-disigned-f or  temperatures  can  arrise  in 
these  regions  and  may  impare  the  safe  operation  of  tie  facility.  These 
regions  of  relatively  high  temperatures  arrise  from  two  principle  causes: 

(i)  because  of  the  nonuniformity  of  the  porosity  distribution, £ ,  regions 
of  relative  minima  in  £  develope  where  relatively  large  heat  transfer  rates 
occurr  between  solid  and  gas,  (thereby  resulting  in  higher  temperatures), 
and  (ii)  because  of  blockages  that  may  develope  as  a  result  of  a)  the 
retention  of  dust  or  foreign  matter  or  b)  of  the  breakage  of  particles  due 
to  possible  t vermal  stresses,  higher  temperatures  may  result  as  well. 

It  would  therefore  be  useful  to  be  able  to  mathematically  model  the 
problem  of  heat  transfer  in  a  bed  of  active  particles  where  variations  in 
the  porosity  occurr  or  blockages  exist,  and  subsequently  develope  a  numerical 
technique  for  distributions.  Such  is  the  purpose  of  this  work. 

The  paper  will  be  divided  into  six  chapters:  the  mathematical  state¬ 
ment  of  heat  exchange  in  an  active  packed  bed;  numerical  approach  used  for 
the  solution;  optimization  and  verification  of  numerical  approach;  study 
of  the  formation  of  hot  spots  due  to  variations  in  the  porosity  distri¬ 
bution;  study  of  formation  of  hot  spots  due  to  blockages;  and  summary  of 
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Important  results  and  conclusions. 


1.  Mathematical  Statement  of  the  Problem 

The  first  choice  that  must  be  made  in  modeling  heat  transfer  in  a 
packed  bed  is  whether  to  enforce  the  conservation  of  energy  by  considering 
the  gas/solid  mixture  or  by  considering  the  gas  and  solid  separately;  i.e. 
to  consider  a  one  phase  or  two  phase  system.  A  two  phase  approach  is 
chosen  for  this  work.  In  addition,  as  variations  in  the  ft  direction  are 
assumed  to  be  zero,  an  axis-symmetric,  two-dimensional  system  is  used. 

The  pair  of  differential  equations,  conserving  energy  in  the  gas 
and  solid  state  respectively,  are  derived  in  appendix  A.  The  results  are: 

conservation  of  energy  in  the  gas: 

(I.Up'.Cp)^  +  i  )  .  i.i  ,  .  hy •  (Ts  -  T 

(?.'  conservation  of  energy  in  the  solid: 

(2.)  -  ^  ^  (r-ks*.*^  )  ■  hy-CTg  -  Ts)  +  Q'*(!  -fc) 

The  basic  assumptions  are  enumerated  in  appendix  A.  ks*  and  kg*  are 
effective  thermal  conductivities.  It  is  left  to  express  these  quantities 
in  terms  of  their  molecular  counterparts  ks  and  kg.  kg*  is  given  as  in 
reference  (  1  ):  kg*  =  £-kg  +  This  form  is  identical,  save  for  th 

factor  £,  to  what  the  effective  thermal  conductivity  becomes  when  one 
assumes  an  eddy  diffusivity,  associated  with  the  Reynolds  stresses  in  the 
averaged  turbulent  momentum  equation,  that  is  defined  analogously  as 
is  defined  for  the  viscous  stresses  in  the  Navier-Stokes  equations. 

kg*  has  an  additional  term,  p.Cp*E,  which  takes 
account  of  the  heat  dissipated  by  the  turbulent  f lucuations  .  We  take  ks* 
as:  ks*  "  (1  -&)*ks.  However,  this  form  for  ks*  only  takes  account  of 
the  gross  effect  that  porosity  may  have  on  the  effective  thermal  con¬ 
ductivity.  It  fails  to  consider  that  heat  conducted  through  the  solid 
phase  really  depends  upon  the  number  of  contact  points  between  spheres  in 
a  given  region.  Given  that  conductivity  based  on  the  number  of  contact 
points  is  hard  to  define,  we  use  the  above-mentioned  expression  for  ks*. 

It  is  important  to  notice  that  theae  expressions  for  kg*  and  ks*  take 
account  of  the  fact  that  not  all  of  the  surface  of  any  given  face  of  V 
(the  differential  volume  in  appendix  A.),  is  entirely  available  to  either 
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the  conductive  energy  flux  of  the  gas  phase  or  of  the  solid  phase.  This 
is  expressed  by  the  terms  6/kg  and  (1  -£)-ks. 

In  the  construction  of  the  differential  equqtions,  we  make  the 
assumption  that  the  porosity  "at  a  point",  based  on  an  infinitessimal 
volume  element,  is  the  same  as  what  would  be  given  if  the  porosity  were 
based  on  any  planar  side  of  the  infinitessimal  volume  element;  i.e.,  at 
a  given  point,  phased  on  volume,  is  equal  to  £,  based  on  area.  However 
talking  about  porosity  defined  at  a  point  is  rather  misleading;  for, 
given  a  point  at  r,©,x,  strictly  speaking  the  porosity  is  either  zero  or 
one.  Therefore  the  porosity  is  discontinuous  and  the  application  of 
partial  differential  equations  is  invalid.  But,  as  our  partial  differ¬ 
ential  equations  do  not  incorporate  variations  in  the  ©  direction,  we 
may,  hypothetically  anyway,  average  out  the  variations  of  £  in  the  © 
direction.  Assuming  such  a  function  could  be  experimentally  derived, 
this  yields  a  continuous  function  in  x  and  r.  See  figure  2.  for  a 
clarifying  illustration.  A  function  in  r  that  we  use  in  much  of  our 
calculations  is  as  given  in  reference  (3)5 

(3.)fc=  0.4  +  0.3  Exp  (-3.5  (g  (1  -  ~>)  +  cos  (6.28  (^  (1  -  |j))> 

Exp  (-0.8  (g  (1-|£»> 

Equations  (1.)  and  (2.)  form  a  coupled  elliptic  system  for  the  un¬ 
known  temperatures  Tg  and  Ts .  Therefore  Diriclet,  Neumann  and/or-mixed 
boundary  conditions  must  be  specified  for  both  the  solid  and  gaseous 
temperatures,  completely  around  the  boundary.  The  boundaty  conditions 


used  are: 


(x  =  0,  r) :  Tg  =>  To; 


0,  or  Ts  =  To 


(x,  r  =  RT):  -kg*  =  hgw»(Tg  -  Ta)  ;  -ks* 
(x  =*  HT,  r):  -  0 

0X  0X 


hSw* (Ts  -  Ta) 


In  passing,  the  other  parameters  used  in  (1.)  and  (2.)  include  hy»  the 
volumetric  heat  transfer  coefficient  which  governs  the  rate  of  heat  trans¬ 
fer  between  solid  and  gas,  and  Q'  the  rate  of  heat  production  per  unit 
solid  volume. 

The  boundary  conditions  at  the  inlet  reflect  (i)  that  the  gas  enters 
at  a  known  temperature  and  (ii)  that  the  solid  could  also  maintain  this 
temperature  or  could  assume,  as  the  gas  and  solid  temperature  fields  do 


at  the  exit,  what  we  believe  the  more  realistic  condition  that  the  normal 
derivative  of  the  temperature  is  zero.  The  boundary  conditions  at  the  long¬ 
itudinal  axis  simply  reflect  an  assumed  symmetry  property.  The  boundary 
condition  at  the  wall  is  a  mixed  one.  Setting  h^^/kg*  and  hsw/ks*  equal 
to  zero  or  a  very  large  number,  an  adiabatic  or  an  isothermal  wall  condition, 
respectively,  is  enforced. 

2.  Numerical  Treatment  of  the  Problem 

In  preparing  (1.)  ,  42.)  and  (4.)  for  numerical  treatment,  it  is 
conenient  to  non-diraensionalize  them  beforehand.  Using  the  expressions: 


x*  »  x/HT,  r*  -  r/RT,  U*  **  U/Vo,  V*  =>  V/Vo,  T  =  Tg/To,  8  -  Ts/To, 
Pr  =  JJc Cp/kg,  Re  =*  Vo*RT/V  ,  Nu  =■  h-Dp/kg,  ks*  «*  (1  -£)*ks, 
kg*  =  £*kg  +  ^*Cp*E 

we  get  the  following  non-dimensional  equations  and  boundary  conditions 
-1 


(5.) 


(6.) 


(7.) 


_  RT 

Re  *Pr  HT2"  ^  fcx*2 

,U*  -1  .RT 

n2 


(fc+Pr-(E/v»^.2  + 


(5 


HT  -  Re*Pr  HT 
.VI  -1  RT 
lRT  "  Re-Pr *RT 
6;(1  -fi,)  Nu.RT  . 

Re.Pr  Dp2  (®  T) 

l2 


■  +  + 

.  p_.A  /F/t,n  .  I  (fc+  Pr-(E/y))-RT  frT 
2  W  +  W>  r*  Re* Pr* RT^  >  *>r* 


;RT  .  . 


HT 


Re*  Pr 


(kS/kg)^f,2  :  a,.il_^l.(ksAs)^2)i2  + 


RT 


RT  • Re  -Pr 

6-ji 


.(( 


& 


(ks/kg)  -  (1 


Nu* 


RT 


Dp^  Re .  Pr 

(x*  =*  0,  r*) :  T  =  1, 


(T 


<x\  r* 


1,1 


lx*  =  1  r*)  • 
ix  -  X.  r  ).  dx» 


(x*,  r*  *  0) 


.  dl 

'  dr* 


e)  -,RT-Q':X1  -fr> 

Re-Pr* kg. To 
K  $ 

“  0  or  0  a  1 

Jx* 

(hgw-RT/kg*).(T  -  Ta/To) 
(hsw’RT/ks*) • (8  -  Ta/To) 
0 

0 


dx* 
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Such  a  system  of  equations  and  boundary  conditions  must  in  general  be 
solved  numerically.  A  constant  space  grid  that  is  used  is  indexed  below: 
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Vriting  (5.)  and  (6.)  symbolically  as 

-  “1^.2  +  U1  +  V1  +  H'T  ‘  "•« 

"M2  »f*2  "  002  if*2  +  U2  |f*  +  v2  |f*  +  "•«  ‘  »’T  +  Q 

and  given  the  finite  difference  expressions, 

^J*2>ij  =  (Tij+1  -  2  Tij  +  Ttj-l)/DR2 

^”*2>ij  3  (Ti+lj  -  2  Tij  +Ti_lj)/DX2 
(|j*>ij  -  (Tij  -  Ti_ij)/DX 

)  i  j  =  (Tij+i  -  Tijii)/2*DR,  simimarily  for  8,  one  can  derive 
the  following  difference  equations  for  T  and  8: 

<’•>  TiJ-i  <5F  -  +  ^  +  +  S  +  H)  +  CW  +  hj 

etj  (H)  -  Ti^j  <££  +  §>  +  Ti+1J  (*p> 

,n  1  n  ,-CC2  V2 .  ,  „  .  ,AA2  CC2.  ,  t'2  ,  .  ,  „  ,-CC2  ,  V2 _ . 

(  .)  ®ij-l  (  ^7  -  2.DR)  8ij  (2(^x?  +  dr7  +  DX  +  ^  eij+l  (  dr2"  +  2-DR^ 

_  a  ,AA2  .  l.'2v  .  _  ,AA2v  ,  _ 

-  lij  (H)  »  t->i-lj  (^77  +  ^  +  *i+lj  (H?>  +  Qij 

One  could  consider,  in  the  solution  of  (8.)  and  (9.),  that  temperatures 
at  all  nodal  points  are  unknown  values .  But  as  this  would  lead  to  the 
inversion  of  an  extremely  large,  sparse  matrix,  an  alternate  method 
should  be  found.  The  problem  in  employing  any  alternate  scheme  is  the 
existence  of  the  terms  in  (8.)  and  (9.)  contributed  by  the  axial  con¬ 
duction;  namely  Ti+jj  and  6i+ij  .  If  the  temperatures  along  row  i  are 

to  be  solved  for,  then  guessed  values  must  be  assumed  for  these  terms. 

An  iterative  procedure  is  clearly  called  for.  This  is  the  complication 
introduced  by  the  elliptic  nature  of  the  equations. 

Assuming  t hat  all  temperatures  are  known  in  tie  kill  iteration  at 
all  grid  points  up  to  and  including  the  (i-1)— -  row  and  in  the  (k-1)  — 
iteration  in  t  ie  iEll  row  and  above,  t  e  coupled  differenced  equations 
(g.)  and  (9.)  are  applied  to  tie  temperatures  in  the  i£ii  row  of  grid 
points  for  all  j .  A  special  type  of  matrix  equation  is  obtained  whose 
form  is  s  iown  in  figure  3,  Employing  an  efficient  subroutine  to  solve 
this  matrix  system  for  T^j  and  Bjj  along  row  i,  t  ese  values  are  sub¬ 
sequently  "over-relaxed"  by  the  use  of  a  met  od  which  employs  the 


T 
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following  formulas: 


(Tij)k  =  (Tipk-1  +  ’’S  *((Tij)matrix  -  (Tlj)k-1^ 

C«i j)k  =  (*ij>k-l  +  ,JS*((«ij)matrlx  -  (Sij)k-l) 


T  lis  procedure  Is  termed  relaxation  by  lines  and  Is  used  to  accelerate 
t  e  convergence  of  the  solution  method.  WS  is  called  the  over -relaxat ion 
factor,  ©hose  value  is  typically  between  1  and  2.  As  we  proceed  to  the 
next  row,  we  can  use  the  newly  updated  values  of  the  i£j.  row  in  our 
calculation  of  t  le  values  in  the  (i+l)st  row  as  given  by  the  matrix  eq¬ 
uation.  Proceeding  in  this  manne”  through  the  bed  one  can  check  whether 
the  largest  error  between  iterations  satisfies  some  error  criterion.  If 
not,  one  must  start  the  process  again,  repeating  as  many  iterations  as 
necessary  so  as  to  obtain  convergence. 

As  a  final  note,  it  may  have  been  observed  that  the  finite  difference 

expressions  used  for  v  in  deriving  (I.)  was  a  first  order,  one-sided 

Ox  .  _ 


difference  expression;  unlike  the  central  difference  eauation  for 
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The  reasoning  for  t:,is  lies  in  t '.e  assertion  that  for  high  enou  g  •  He  tie 
first  derivative  term  in  equations  (5.)  is  significantly  greater  in 
magnitude  than  other  terms,  and,  according  to  reference  (4),  to  avoid 
numerical  instabilities  in  the  solution  as  well  as  enhance  the  accurracy 
of  the  solution,  the  one-sided,  or  upwind,  differencing  scheme  should  be 
employed  to  .  Unfortunately  the  upwind  differencing  scheme  was  also 
applied  to  C”  .  As  no  relative  weighting  with  respect  to  Re  exists  be¬ 
tween  the  first  and  second  derivative  terms  in  equation  (2,),  the  rational 
for  using  the  upwind  schemefgf  was  inappropriate  for  use  on  ^  . 
Therefore  it  was  inconsistent  to  apply  upwind  differencing  to  and 
it  is  possible  that  the  solutions  obtained  were  not  as  accurrate  as  could 

have  been  obtained  using  central  differencing  on  ;  however  the 

0'-< 

discrepencywould  be  small  and  would  have  no  implications  on  conclusions 
to  be  drawn. 


3,  Optimization  and  Verification  of  Numerical  Method 

This  chapter  concerns  itself  with  the  actual  iraplimentation  of  the 
numerical  method  in  obtaining  solutions  with  the  purpose  of  minimizing 
the  number  of  iterations  necessary  for  convergence  and  the  comparison  of 
the  actual  results  with  those  generated  bya  priviously  established  num¬ 
erical  proceedure  which  Ignores  axial  conduction.  Besides  the  usual 
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geometric  and  physical /thermal  properties  of  the  gas  and  solid,  the  porosity, 
the  velocity  field,  the  Nusselt  number  and  eddy  diffusivity,  E,  must  be 
specified  in  order  to  generate  a  solution,  t’sing  a  reasonable  value  for 
Nu;  setting  E  to  zero;  and  as  required  by  the  analytical  solution  , 
assuming  the  velocity  and  porosity  are  constants,  (V/Vo  -  1  and  =  0.4), 
the  numerical  proceedure  as  presented  in  this  paper  was  applied.  As  for 
boundary  condition,  the  gas  and  solid  were  assumed  at  a  given  temperature 
at  a  given  temperature  at  the  inlet  as  well  as  at  the  wall. 

Although  the  discretization  method  as  described  by  equations  (8.)  and 
(9.)  was  found  to  be  quite  satisfactory,  in  our  first  attempts  to  achieve 
a  solution,  an  alternate  method  was  also  considered.  This  alternate 
method  is  equivalent  to  employing  (8.)  and  (9.)  except  that  the  nodal 
temperatures  multiplying  the  coefficients  Ul,  l’"’,  Vl,  V?,  were  considered 
at  the  (k-l)st  iteration.  The  obvious  superiority  of  the  method  represented 
by  O.)  and  (9.)  is  reflected  in  figure  4  .  The  advantages  are  shown  in 
two  ways:  (i)  a  much  lower  number  of  iterations  is  necessary  for  con¬ 
vergence  and  (ii)  the  first  method  is  less  sensitive  to  changes  in  WS  than 
the  second.  what  may  also  be  seen  in  figure  4.  is  an  agreement  with  the 
assertion  that  for  an  elliptic  equation  with  constant  coefficients,  the 
WS  that  results  in  the  minimum  number  of  iterations,  lies  between  1  and  ?. 
This  is  seen  to  hold  even  for  a  variety  of  Re. 

The  method  of  calculation  having  been  established,  it  remains,  as  a 
further  check,  to  compare  the  results  of  the  numerical  solution  presented 
here,  (assuming  constant  porosity,  velocity  and  an  isothermal  wall  con- 
cition),  to  those  results  given  by  an  associated  numerical  solution 
which  has  already  been  developed  and  which  ignore*  axial  conduction. 

This  numerical  proceedure  is  based  on  the  parabolic  mathematical  equations 
which  one  obtains  when  axial  conduction  is  ignored.  Such  a  system  has 
the  numerical  advantage  that  the  solution  can  be  found  directly  without 
Iteration.  The  results  of  the  numerical  solutions  for  a  Re  of  five  thousand 
are  shown  ip  figure  5.  In  addition,  the  results  of  the  previous  num¬ 
erical  solution  ignoring  axial  conduction;  are  compared  to  an  analytical 
solution  which  ignores  axial  conduction  by  construction  as  well  and  which 
can  be  generated  in  terms  of  Bessel  functions.  This  comparison  is  almost 
exact.  As  the  difference  between  the  numerical  proceedures  is  the  in¬ 
clusion  of  axisi  conduction,  it  follows  that  the  difference  in  the  num- 


erical  results  is  due  to  the  conduction  of  heat  in  opposition  to  the  dir¬ 
ection  of  the  flow.  Looking  at  the  non-dimensional  equation  (5.),  one 
sees  that  the  product  Re*Pr  determines  the  relative  importance  of  the  con¬ 
vective  terms  relative  to  the  conductive  terms .  Therefore  as  the  Re  in¬ 
creases,  one  would  expect  that  the  diffusive  terms  play  an  increasingly 
less  significant  role  and  as  figure  6.  shows,  that  the  numerical  solutions 
should  approach  one  another.  Thus  for  Re*Pr  down  to  about  "50  the  solution 
~iven  by  the  numerical  proceedure  including  axial  conduction  should  approx¬ 
imate  the  numerical  solution  excluding  axial  conduction  to  within  17,  . 

This  implies  that  there  exists,  for  given  parameters,  a  Re,  above  which  the 
simpler  parabolic  system  of  equations  can  be  used;  and  below  which  the 
more  complex  elliptic  system  must  be  employed.  However,  as  seen  in  figure  7. 
this  relative  error  does  not  exhibit  a  uniform  distribution,  but  is  a 
function  of  the  position  of  the  bed.  The  result  that  the  relative  maxima 
of  error  occurr  near  the  exit  and  near  the  wall  of  the  bed,  can  be  ex¬ 
plained  by  figure  8.  .  Figure  8.  shows  that  the  second  derivative  of  the 
gas  temperature  is  highest  near  the  bed  exit  and  near  the  bed  wall;  and 
therefore  one  would  expect  that  in  these  regions  the  relative  discrepency 
between  the  numerical  solutions  would  be  greatest . 

It  is  seen,  therefore,  that  the  numerical  proceedure  described  in  this 
paper  can  be  optimized  in  the  sense  that  the  number  of  iterations  are  min¬ 
imized,  by  keeping  WS  between  1  and  1.2  .  In  addition,  the  agreement 
between  the  results  of  the  numerical  approach  presented  here  and  the  results 
of  the  numerical  proceedure  which  ignored  axial  conduction,  being  shown 
to  become  better  as  the  Re  increases,  serves  to  justify  that  the  program 
is  operating  properly. 

4.  Study  of  the  Formation  of  Hot  Spots  Due  to 
Variations  in  the  Porosity  Distribution 

In  order  to  better  simulate  the  local  properties  of  a  packed  bed 
and  its  heat  transfer  characteristics,  experimental  expressions  for  the 
porosity,  given  as  equation  (3.),  for  the  local  Nusselt  number,  given  in 
appendix  Cl,  and  expressions  for  parabolic  and  turbulent  superficial  vel¬ 
ocity  profiles  are  used.  The  expressions  used  for  the  velocity  profiles 
are  given  in  appendix  C2  and  are  graphed  along  with  the  porosity  and 
Nusselt  number  distributions  in  figures  9.  and  10.  respectively.  The 
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oscillatory  variation  of  the  porosity  is  a  result  of  the  influence  of  the 
wall  containing  the  packed  bed.  The  porosity  is  1  by  definition  at  the 
wall  and  as  the  center  is  approached  the  porosity  becomes  a  constant 
value,  (typically  0.4  for  a  lightly  packed  random  bed).  Figure  11.  shows 
the  gas  temperature  obtained  at  the  exit  of  the  bed  for  a  constant  velocity 
profile  and  type  (i)  boundary  conditions,  (see  '>clpv  )  .  The  oscillatory 
variations  observed  are  due  to  the  variations  in  the  porosity  distribution. 
This  can  be  explained  by  realizing  that  regions  of  relatively  low  porosity 
are  regions  of  high  spherical  density  and  therefore  regions  where  high 
heat  transfer  rates  between  the  solid  and  the  gas  exist.  The  dashed  line 
in  figure  11.  ,  is  the  numerical  solution  given  by  the  program  that  ignores 
axial  conduction.  It  shows  good  agreement  with  the  present  solution. 

Figures  12.  and  13.  show  the  solution  obtained  at  the  bed  exit  for  assumed 
turbulent  and  laminar  velocity  profiles  respectively.  Type  (i)  boundary 
conditions  are  used  here  also.  Again  the  oscillatory  variations  are  due 
to  the  local  variation  in  porosity;  but,  in  addition,  due  to  the  fact  that 
the  velocity  is  greater  than  the  mean,  up  to  a  certain  distance  from  the 
center,  and  is  less  than  the  mean  nearer  the  wall,  the  cooling  affects  of 
the  fluid  are  more  efficient  nearer  the  center  than  nearer  the  wall.  There¬ 
fore  the  relative  maxima  in  the  temperature  profiles  are  due  to  both  the 
relative  minima  of  porosity  that  occurr  near  the  wall  and  to  the  inefficient 
cooling  in  this  region.  It  should  be  mentioned  that  E  was  assumed  zero  in 
these  cases  involving  figures  11.,  12.,  13.  . 

Assuming  a  turbulent  velocity  profile,  figure  14.  presents  the  re¬ 
sulting  gas  temperature  distributions  for  a  variety  of  boundary  conditions . 
The  boundary  conditions  used  are  defined  as  follows: 

type  ’  (i):  Ts(0,r)  ■  To,  Ts  -  Tg(x,RT)  »  Ta 

type  (ii):  (0,r)  *  0,  Ts  ■  Tg(x,RT)  ■  Ta 

type  (iii) :  ^  (0,r)  "  0,  ^  (x,RT)  -  0  ; 

other  boundary  conditions  remain  as  previously  described 
and  are  the  same  for  each  of  these  cases. 

It  can  be  seen  that  by  changing  from  a  type  (i)  to  a  type  (ii)  boundary 
condition,  the  temperature  profile  of  the  gas  increases  yet  keeps  the 
general  shape  of  the  isothermal  case  and  that  the  position  of  the  hot 
spot  has  kept  the  same  position,  with  the  introduction  of  the  type  (ill) 


boundary  condition,  the  opportunity  for  heat  to  be  passed  through  the  wall 
is  lost.  Therefore,  as  expected,  the  temperature  profile  is  raised  even 
further.  In  addition,  the  hot  spot  is  moved  significantly  closer  to  the 
wall . 

It  should  be  mentioned  in  passing,  that  as  with  the  results  in  figures 
11,,  12.,  and  13.  the  dashed  line  in  figure  14.,  is  calculated  using  a  Nusselt 
number  bqsed  on  the  mean  Re;  whereas  the  other  curves  in  figure  14.  are 
calculated  using  a  Nusselt  number  based  on  the  local  Re.  The  shift  to  a 
local  Nu  slightly  raises  the  temperature  profile  due  to  the  fact  that  the 
local  Re  (and  therefore  the  local  Nu),  is  greater  than  the  mean  Re, 

(mean  Nu) ,  for  most  of  the  radial  dimension.  This  allows  a  larger  heat 
transfer  between  the  solid  and  the  gas  which  raises  the  resulting  temper¬ 
ature  profile. 

Earlier,  the  results  of  the  numerical  technique  developed  in  this  paper 
compared  to  a  previously  established  numerical  technique,  (ignoring  axial 
conduction),  which  in  turn  was  seen  to  be  identical  to  an  associated  ana¬ 
lytical  solution.  The  results  of  the  numerical  proceedures  were  seen  to 
approach  one  another  as  the  Re  increased.  This  implied  that  the  numerical 
proceedure  including  axial  conduction  was  implimented  properly .  To  attempt 
to  justify  the  numerical  results,  one  must  compare  them  with  experimental 
ones.  Such  experimental  work  was  done  at  the  Von  Karman  Institute  and 
is  discussed  in  reference  ( 5 ) .  Using  their  given  data,  numerical  cal¬ 
culations  were  made.  The  experimental  and  numerical  results  appear  in 
figure  15 .  .  As  the  temperatures  are  scaled  with  respect  to  the  hot  spot 
temperature,  only  the  predicition  of  the  position  of  the  hot  spot  can  be 
made  and  not  its  magnitude.  Considering  that  the  experimental  data  was 
in  the  form  of  a  temperature  map  of  the  bed  exit  and  that  one  necessairly 
must  choose  the  radial  profile  to  compare  his  numerical  results  with, 

(the  radial  profile  through  the  hot  spot  of  the  bed  exit  was  taken) 
and  that  one  cannot  expect  the  porosity  of  the  experimental  packed  bed 
to  be  the  same  as  that  used  in  the  numerical  program,  that  the  numerical 
results  rather  successfully  predicted  the  radial  position  of  the  hot  spot 
as  well  as  the  profiles'  general  shape,  cap  be  inferred.  An  additional 
point  should  be  mentioned  concerning  figure  15,  .  As  the  bed  height  is 
Increased,  the  level  of  temperature  at  the  exit  is  raised  and  therefore 
the. greater  the  heat  loss  through  the  wall.  This  accounts  for  the  inward 
movement  of  the  hot  spot . 


Tt  is  seen,  then,  that  the  variation  in  porosity  introduces  hot  spot 
regions  that  exist  close  to  the  wall  of  the  packed  bed  and  that  assumed 
parabolic  or  turbulent  velocity  profiles  amplify  the  magnitude  of  these 
hot  spots.  In  addition,  the  affect  of  employing  adiabatic  wall  conditions 
tends  to  increase  even  more  the  temperature  of  these  regions. 

As  a  final  point,  the  results  presented  in  fighures  14.  and  15.; 
unlike  figures  11.,  12.,  13.,  were  calculated  using  a  value  of  20  for  the 
modified  Peclet  number,  Pe'.  This  value,  dependent  upon  the  Re  and  the 
ratio  Dp/RT,  is  obtained  from  the  experimental  data  given  in  reference  (  i  ) . 
Because,  (i)  the  value  of  Pc*  determines  the  etddy  diffusivity,  E,  which 
in  turn  contributes  to  the  effective  thermal  conductivity  kg*,  and  (ii) 
the  greater  the  value  of  kg*,  the  greater  the  dissipative  capacity  of  the 
gas  is  to  damp  out  rapid  temperature  variations,  it  is  important  to  be 
able  to  represent  Pe'  satisfactorily  if  one  wants  to  predict  the  magnitude 
of  hot  spots.  As  an  example,  the  calculated  curve  of  figure  15.,  for 
HT/Dp  =  10,  is  repeated  in  figure  16.  and  compared  to  a  solution  calculated 
from  the  same  data  save  that  the  Pe'  is  very  large,  (therefore  implying 
E  =*  0)  . 

As  the  assumed  axial  velocity  profile,  or  for  that  matter  the  "true" 
local  mean  axial  velocity  profile,  would  not  be  expected  to  show  great 
variations  in  the  radial  direction,  which  would  cause  relatively  large 
eddy  diffusivities  to  arrise,  the  use  of  a  constant  E  seems  reasonable. 
However  if  one  assumes  a  blockage  exists,  large  perturbations  in  the  flow 
field  would  result;  thereby  causing  large  gradients  and  therefore  larger 
eddy  diffusivities.  In  addition,  as  there  would  be  little  mass  flow  in  the 
vicinity  of  the  blockage,  the  conductive  terms  would  become  the  main  heat 
transfer  process,  whatever  the  bulk  Re.  Therefore  only  the  application  of 
the  numerical  proceedure  including  axial  conduction,  as  presented  here, 
would  be  justified. 


5.  Study  of  the  Fornaticn  of  'lot  Spots 
Duo  to  the  Presence  of  blockages 

In  the  operation  of  a  packed  bed  heat  exchanger,  a  blockage  concicvnhly 
might  result  from  the  retention  of  dust/foreign  matter  being  convrcted  by 
the  gas.  In  the  region  of  the  blockage,  both  the  inefficient  cooling  by 
the  gas  and  the  heating  of  the  dust  surrounding  the  active  spheres,  would 
lead  to  the  formation  of  a  hot  spot.  The  blockage,  itself,  is  introduced 
into  the  program  as  a  region  of  the  relatively  low  porosity  of  0.1  .  Out¬ 
side  the  blocked  region  the  porosity  is  0.4  .  To  avoid  large  gradients, 
a  linear  variation  in  the  porosity  between  the  blockage  and  the  rest  of  the 
bed  is  used.  See  figure  16.  for  a  clarifying  illustration.  The  local  vel¬ 
ocity  vector,  (axial  and  radial  components),  is  calculated  by  a  program 
that  solves  the  non-dimensional  form  of  the  Ergun  equation.  Refer  to  ref¬ 
erence  (6)  for  the  equation  and  examples  if  its  applications.  The  velocity 
field,  for  the  case  of  the  blockage,  is  shown  in  figure  17.  .  As  expected 
the  flow  tends  to  be  diverted  around  the  blockage,  thereby  leaving  the 
blocked  region  with  little  mass  flow  through  it . 

Unlike  the  previous  cases  where  E  was  assumed  constant  over  the  entire 
field,  we  would  expect  relatively  large  velocity  gradients  in  the  case  of 
the  blockage  and  therefore  regions  of  higher  eddy  diffusivities .  A  velocity 
gradient  model  relating  E  to  the  velocity  field  is  used  and  is  explained 
in  appendix  D. 

Using  the  velocity  gradient  model  for  eddy  diffusivity,  the  above 
definition  for  porosity,  the  Nusselt  number  expression  in  appendix  Cl  and 
the  calculated  velocity  field,  temperature  distributions  were  obtained 
for  a  variety  of  cases . 

Figure  IT.  shows  axial  and  radial  gaseous  temperature  curves  for 
constant  as  well  as  variable  E,  using  an  isothermal  wall  condition.  It 
must  be  mentioned  that  in  calculating  the  variable  eddy  diffusivity  dis¬ 
tribution  from  the  velocity  field,  an  arbitrary  constant  multiplying  the 
vector-  components  of  E,  Ex  and  Er,  was  adjusted  so  that  the  vectoral  sum 
of  the  average  Ex  and  the  average  Er  over  the  entire  bed  mesh,  was  such 
that  the  effective  Pe’  calculated  from  this  mean  value,  matched  the  Pe' 
given  in  reference  (  .  As  figure  IS.  shows,  the  presence  of  a  blockage 

causes  a  signigicant  temperature  peak  in  the  radial  profile  and  effects 
an  acceleration  in  the,  what  would  otherwise  be  a  linear  temperature  profile. 


-13- 


The  variable  E  solution  save  larger  estimates  of  peak  temperatures  than 
the  constant  E  solution.  This  can  be  explained  by  noting  in  figure  19,  that 
in  the  region  of  the  peak  of  the  radial  temperature  curve,  the  E  estimated 
by  the  variable  model  is  much  lower  than  the  mean  bed  value.  Hence  less 
turbulence  and  less  efficient  cooling  is  the  result  and  therefore  the  var¬ 
iable  E  model  will  give  greater  temperatures.  Identical  arguments  can  ex¬ 
plain  the  discrepency  between  the  axial  temperature  profiles  of  figure  IS. 
when  one  notes  that  as  with  the  upper  curve  in  figure  19.,  the  lower  one 
shows  that  the  variable  E  is  much  less  than  the  mean  value  in  the  region 
of  interest . 

As  mentioned  in  chapter  ,  the  expression  ks*  =  (1  -£,)*ks,  does  not 
take  into  account  that  ks*  really  depends  on  the  number  of  contact  points 
between  spheres  rather  than  exclusively  on  the  porosity.  One  would  expect, 
then,  that  ks*  would  in  reality  be  somewhat  less  than  that  given  by  the 
above  expression.  Figure  20.  shows  the  radial  temperature  at  the  bed  exit 
for  the  same  blockage  case  as  treated  before,  (given  as  dots  (•)),  and 
compares  this  result  with  those  profiles  calculated  by  assumming  subse¬ 
quently  lower  values  of  ks.  Because  of  increasing  resistence  to  conductive 
heat  flow,  the  model  predicts  subsequently  larger  peak  temperature  values 
as  ks  decreases. 

The  affect  of  employing  a  type  (iii),  adiabatic,  wall  condition  as 
opposed  to  an  isothermal  one,  has  virtually  no  affect  on  the  character¬ 
istics  of  the  mid-bed  radial  profile,  as  exemplified  by  the  upper  curve 
in  figure  21.  .  Not  only  is  this  ineffectuality  in  changing  pfeak  temper¬ 
ature  values  true  at  the  mid-bed  position,  but  it  is  true  for  the  entire  bed. 

The  lower  curve  in  figure  21.  shows  that  as  the  blockage  is  moved 
closer  to  the  bed  exit,  the  temperature  peak  at  the  bed  rises  accordingly. 
This  brings  up  the  immediate  concern  that  temperature  peaks  may  develope 
inside  a  packed  bed  with  a  blockage  that  may  exceed  the  temperatures 
realized  at  the  bed  exit;  and  as  it  is  difficult  to  determine  temperatures 
inside  an  operating  packed  bed  heat  exchanger,  these  temperatures  may  go 
undected.  The  results  from  which  the  lower  curve  in  figure  21.  was  derived. 
Justify  this  concern;  for  as  the  blockage  is  in  either  the  lower,  middle, 
or  upper  position,  the  maximum  gaseous  temperature  in  the  blocked  region 
attains  39.5%,  (at  X  =»  0.47  HT?,  95 .5%,  (at  x  =  0.63  HT)  and  99.9%,  (at 
x  »  C.°0  HT),  of  the  maximum  value  of  the  exit  temperature  respectively. 


r/ 


m;  / 


i 
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It  is  therefore  seen,  then,  that  the  introduction  of  blockages  into 
packed  beds,  result  in  large  peaks  in  the  temperature  field  distribution. 

Sy  using  a  velocity  gradient  model  for  the  eddy  diffusivity,  E,  the  temper¬ 
ature  peaks  have  br.en  higher  than  might  have  otherwise  been  expected  had  a 
constant  E  been  used.  These  hot  spots  may  attain  higher  values,  still, 
if  one  assumes  that  the  effective  thermal  conductivity  of  the  solid  is  lower 
than  that  given  by  the  assumed  expression.  The  magnitude  of  these  peaks, 
in  addition,  are  functions  of  the  position  of  the  blocked  region;  and  it 
was  seen  that  the  possibility  definitely  exists  that  the  magnitude  of  hot 
spots  may  be  greatest  inside  the  bed,  itself,  and  therefore  difficult 
to  detect. 

6.  Summary  and  Conclusions 

Envisioning  the  use  of  an  active  packed  bed  in  a  nuclear  reactor, 
the  concern  of  hot  spot  development  due  to  porosity  variation  or  due  to 
the  presence  of  blockages,  is  raised.  It  would  be  useful,  therefore,  to 
be  able  to  calculate  the  resulting  temperature  field  assuming  that  these 
porosity  variations  or  blockages  occurr.  To  that  end,  a  two  phase  axis- 
symmetric  mathematical  model  for  heat  transfer  in  an  active  packed  bed 
was  proposed.  As  axial  conduction  is  considered,  the  resulting  equations 
are  elliptic  and  if  calculations  are  to  be  done  efficiently,  these  equations 
must  be  solved  by  an  iteration  proceedure,  A  method  of  over-relaxa'.: .. on  by 
lines  wa s  employed  and  whose  arbitrary  coefficient  of  relaxation  was  chosen 
so  as  to  approximately  minimize  the  number  of  iterations  needed  for  con¬ 
vergence.  A  numerical  proceedure,  which  ignores  axial  conduction,  which 
was  shown  essentially  identical  with  an  associated  analytical  solution, 
gave  results  which  approached  the  numerical  values  of  the  present  program 
as  the  Reynolds  number  increased.  As  the  relative  difference  between 
the  two  solutions  was  on  the  order  of  3-47.  for  low  Re  and  decreased  to 
approximately  17  for  moderate  Re,  the  results  of  the  numerical  solution 
were  verified. 

Concerning  the  case  where  variable  porosity  is  assumed,  it  was  seen 
that  temperature  peaks  developed  close  to  the  bed  wall.  The  position  and 
magnitude  of  these  hot  spots  was  seen  to  be  a  strong  function  of  the  assumed 
velocity  profile,  the  bed  height,  as  well  as  the  boundary  conditions  at  the 
wall.  The  present  solution  was  also  seen  to  be  in  close  agreement  with 
the  results  of  the  previous  program  which  ignores  axial  conduction.  It  is 


in  this  respect  that  Che  following  is  proposed.  For  as  most  industrial 
uses  insure  Re  of  sufficient  magnitude  to  assure  that  the  inclusion  of 
axial  conduction  would  give  temperatures  of  90-95%  of  the  ones  given  by  the 
program  excluding  axial  conduction  and  as  the  elliptic  equations  involve 
more  complicated,  time-consumming  methods,  it  would  perhaps  be  more  efficient 
to  base  a  design  method  on  a  mathematical  model  ignoring  axial  conduction, 
keeping  in  mind  that  the  results  are  to  some  degree  an  upper  limit  to  the 
"actaul  temperatures".  This  concept  of  the  questionable  usefulness  of  the 
inclusion  of  the  axial  conduction,  however,  cannot  be  applied  to  the  case 
where  blockages  are  present. 

As  demonstrated,  even  for  relatively  large  bulk  Re,  the  mass  flow  in¬ 
side  the  blockage  itself  was  small  and  therefore,  in  the  blocked  region, 
heat  transfer  by  conduction  will  play  the  major  role.  The  blockage  intro¬ 
duced  a  substantial  peak  in  the  radial  temperature  profile  through  the 
blockage  itself.  The  overall  maximum  temperature,  however,  was  still  at 
the  bed  exit;  even  though  it  was  approached  by  the  maximum  blockage  temp¬ 
erature  as  the  blockage  moved  downstream.  It  bears  looking  into  the  pos¬ 
sibility  that  the  maximum  blockage  temperature  may  exceed  the  maximum  exit 
temperature  by  effecting  changes  in  the  position  and  size  of  a  blockage. 

As  a  final  point,  it  was  mentioned  that  the  actual  ks*  may  be  lower 
than  assumed.  It  was  demonstrated  that  by  lowering  ks*  the  peak  of  the 
radial  exit  temperature  profile  was  raised.  It  would  therefore  be  useful 
to  investigate  the  true  nature  of  ks*  so  as  one  could  place  more  certainty 
on  whether  a  proposed  design  for  the  operation  of  an  active  packed  bed  is 
acceptable  if  one  assumes  the  presence  of  a  blockage. 
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App  endix  K. 

)  rosraii  Listings  and  Saim.le  Outputs 


Program  code  for  calculating  the  ten-  .  .  . 

perature  field  assuming  variable  velocity 
and  porosity  distributions 

Corresponding  output  . 

Velocity  field  calculation  for  the  case  .  . 

of  a  blockage 

Program  code  for  calculating  the  tea-  .  .  . 

perature  field  for  the  case  of  a 
blockage;  using  the  velocity  d;-ta  above 

Corresponding  output  . 


c 


•  ?6 

.  ?9 

•  3? 

•  4? 

1 7 


r 
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D  E F I HE  FILE  8= 108=N : 10 

C  0  M  '1 0  S'  /  8  L  K  <  /  ft  i  (  4  I  .>  ,  B  1  i  4  1  >  ,  C  I  (  4  1  >  ,  p  1  •;  4  1  )  ,  £  1  <  4  1  ) 

CCMMCNVBL  KE/ AO  (  4  I  ) ,  B  2  I  4  1  )  ,  C  2«  4  1  )  .  0  2 i  4  I  )  .  £  2  <  4  1  ) 

OOMMGN/FZ/~<  43  -43  ),  M  S  U  P  MS  UP 

DIMS  N  S  I  0  S’  V  :  (  ?.  1  -  3  I  )  .  V  2  <  2  i  ,  3  1  )  -  S  CU.  N  1  <  3  I  )  -  S  0  L  N  2  i  3  I  )  ,  £  P  <  3  1  )  -  E  ?  R  •:  3  !  )  , 
I  P  N  <  3  I  >  »  U  3  A  :  3  I  :• 
c 

pe ao  ;  ic5,  son.'  rt ,ht,  re- PP  ,  PNU,  TKC -  TKR 
C  FORMAT*  7F  10  .  0  ) 

v  REAO< 1 05, 20  0  >MMAX»NKAX, J99MAX.EPS, as. DP. TA-OP- POR 
200  f  0RNAT8  31  5- 5F  i  0  0,F7  5) 

READ  •:  105,150)  F  E  -  T  0  ,  V  I  S  ,  M  9  -  R  MG  U  -  R  NS  W  -  R  N  3  0  -  A  L  p  H  A  ,  B  E  T  ft 

15  0  FORMAT  ( 3  <  E 1 0 . 4  > ,  I  3  -  3 <  E 1 0 . 4  ) -  2  < F  4  .  1  >  ) 

r 

•JRITE  I  108,800)  RT, HT . RE . PR , PNU , TKG, TKR , POR, PE » VI  S 
■JR  UE  ■:  10  8,300)  UMAX,  UMAX,  ,!  99  MAX  ,  EPS  ,  WS  ,  DP  ,  TA,  TO,  3P  ,  M9 
3  0  0  FORMAT  <  1 X,  1 0<  El  2 . 4  >, /  ) 

VRITE  < 108,801  >  RNOW, RNSW, RNSO, ALPHA, BETA 
3  01  FORMAT  <  1 X , 5 <  E  1  0  .  4  )  , / /  ) 

C 

C  INITIALIZATION  OF  TEMPERATURES 
DO  5  I  -I,  MM  AX 
DO  5  J  -  1  , N  M  A  X 
'J  1  (  1  ,i  )  ~  ft  L  P  H  A 

a  2  <  i ,  j  >  *  e  e  t  a 

5  CONTINUE 
U  S  =  1  . 

vs=c . 

e  =  dp*r£*vis,-(ft*pe:> 

P  0  R  R  =  0  . 

M  1  =  M  M  A  X  -  1 
N  1  =  N  M  A  X  -  1 
D  X  =  1  .  /  F  L  0  A  T  (  M  1  ) 
f  D  R  =  1  .  /  F  L  0  A  T  (  H  1  ) 

J  9  9  =  0 

S  =  F  0  ft  T  M  i  -  1  >*  F  L  0  AT  <  N  1  -  1  ) 

C  DETERMINATION  OF  COEFFICIENTS  IN  DISCRETIZED  EQUATIONS 
7  J 99=099+1 
DT 1 =0 . 

D  T  2  =  0  . 

DERI  *  0  • 

D  £  R  2  -  0  . 

DO  10  1  =  2 . M 1 
I  P  1  =  I  +  1 
I M  i  =  I  -  l 
DO  20  J  =2 , N 1 
CALL  DATSilT  l  ,  I  SO  ) 

GO  TO  <  827,  829  ),  IS ’J 
829  CONTINUE 

IF  i  I  E3  .  2  >  GO  TO  8  25 
P  0  R  =  E  B  <  J  ) 

PORR=£PR : J ) 

P  N  U  =  P  N  <  0  ) 

u  s  =  u  s  a  <  j  ) 

825  R  =  FLOAT<  J-l  )*0R 

IF  (I  CT .  2  )  GO  TO  8  30 
r 

C  FIND  E ( R ) 

r 

Z  P  =  P  T  * <  1  .  -R  VOP 

POP  =0  4  +  0 . 2 *( EXP ( -2  5  *  Z  R )  +  C  0  S  < 8 . 28*ZR)+£XP< -0  8  *  Z  R  )  ) 

P0pPl-I  05*EXP<-3  5 ♦ Z P  > 

P0RP2- 1 . 884+31 N: 8  28*ZR)*EXP(-0.3fZR) 


I 

c 


c 

r 


C 

c 


t  2  7 


3  f. 


P  2  p y  3  :fj  .  24  *  f  0 St  <  2 '?  * 2  R  j  «  E P '  -  0  3  *  Z  R  ) 
p  o  P  p  =  P  T  «*  ■:  PTRRi  +  P  C  R  k  2  +  P  0  R  R  3  0  /  D  P 
tf!i  > : P 0 2 
CPFS'  2  >  -  r  3  R  R 

R=  FLOAT'  0-  1  MO? 

IF  : M?  . EQ.  I/  50  TO  90 
if  •:  n  ?  e  a  3  :•  go  to  s  o 
y  s  =  2  *  t  1  -  R  *  •*  t  ; 

G  0  TO  3  0 

c  i  r,  -■  q  375  /  ;pr,»3  ?*,  > 

RPL5“?0PT< FLO/ 2 .  > 

REF-PE*P=LP 
y  p  =  r  e  :  *  •:  i  .  -  e  > 

?  ft  -  1  . -30 . /REF 
R B -  1  .-5. /REF 

IF  •:  R  .  G  T  .  R  A  >  GO  TO  S  & 

US--RF_D*7  5  .  5+2 . 5*AL0S<  YP*<  1  .  5*(  1  .  +R  )  )/{  i 
GO  TO  5  o 

IF  <  R  . G  T  .  PS'  GO  TO  ? 7 

y $ - r f :  o * •:  - 3 .  05+5.  *alog<  ypo  > 

GO  TO  30 

y  s  =  R  F .  D  *  r  p 
y  s  ft  •:  ,i )  =  u  3 
reb=re  *d?*us/rt 


FIND  N  U  '  R  E  >  E  I  R  >  > 


+  2  .  *  R  *  R  )  >  ) 


C= 1 . -POR 
0 2  - C  *  *  0 . 3  3  3  3  3 

01-2 . -3 . *02+3. *02**5. -2. *02**6. 

P  N  U  -  i  .  3  *  •'  <  PEB+PR/ IP  0  2*01  )  )  +  (  1  .  -  C  *  *  1  .  6  6  3  6?  )  >  *  *  0  .  33333 
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.  LT  . 

-1  .  39  2E-01 

! 

.  LT  . 

-c*  .  7f  ?E  - 02 

X 

LT  . 

3  .  832E-03 

'  LT.  7.533E-02 

X  . LT .  1 . 4S8E-01 

"  . LT .  2 .  133E-01 

0  . LT  .  2  .  89  8E -  0 1 

.LT.  3.6I3E-01 
?  . LT  .  4  32  SE-01 

,  .LT.  5.043E-01 

M  . LT  .  5  .  758E-01 

!  . LT  .  e  .  47  3E-01 

7  .  I83E-0I 


9 


.  LT  . 


-17- 


0  £  F  I  N  E  FILE  10  =  0  :  1<  R  :  5034  ,  G  :  !  .  BN  ,  F  I  ,  R  U  > 

DEFINE  FILE  5=105=M:SI 
DEFINE  FILE  8M08=M:L0 
COMMGN/DPRH/DP  •  R  T  ,HT 
CO M10S/FIU! PE/DENS. VIS/FE.PR 
C  0  M  M  0  N  /  H  E  ft  T  /  7  K  G  j  T  K  3  ■  Q  P  •  T  A  i  T  C 
COMMON/ Id  ALL  /  RUG  V  >  Q  V  G  ,  R  MS  'J  >  0  US 
COMMON  /  li  T  0  T  /  M  N  ft  >1  ,  H  M  ft  X 

(C  0  M  M  0  N  /  T  C  ft  L  V  /  N  F  ,  N  G  ,  U  S  ,  T  £  S  T 

C  0  N  M  0  N  /  D  ft  T  P  0  R  /'  P  M  1  NftO,  PHAXftO  ,  PM!  N  ft  I  ,  PHftXAl 

COMMON  /  I  ND?  OR  /  11  i  AO  »  I  2  2ft  0  ,  0  t  ’  AC  .  0  4  AC  ,  11  lft  1  ,  I  2?  ft  :  ,  .i  1  1  ft  1  ,  j  4  ft  i 

C  0  M  M  0  S  /  I  N  D  2  I  >  J 

COMMON/ POROUS /POPiPORXiPORP 

C  0  MM  ON  /  I  M  SP/DR  ,  D  X  .•  D  R  2  .  DX  2  ,  DRD  P  >  DXD  X  .  DXDR 

C  0  MM  0 N / B L  K K / A  1 <  4  1  0  .  B  1  <  4  1  )  ,  C  1  <  4  1  )  .  D  I  <  4  1  )  ,  £  1  (  4  1  ) 

CO  MM  ON /PL  KE/ A2<  4  1  T  >  82  <  4 1  >,  C2<  41  )  •  0  2  <  4  1  > , £  2<  4 l  > 

DIMENSION  ‘it  1  (  3  1  ,  3  1  >  .  V  2  <  3  1  .  3  1  >  .  S  OLH  1  <  3  1  >  ,  ?  0  LM  2  !  7  1  1  ,  P  1  3  1  ,  3  1  > 

C 

CALL  LECRIB 
C 

READ  <5, SCO)  OOONftXiEPSiZl  »  P  N  S  0  > ALPHA. -BETA 
BOO  FORMAT  <  15  ,  5<  F 1 0 . 4  >> 

V RITE  <Si?00>  J99NAX, EPS , 2 1  , RNSO  -  ALPHA , BE T A 
7  0  0  FORMAT  <  5  X ,  3  H  J ? 9 M  A  X  =  ,  I  3 , 2  X  ,  5  H  E  P  $ =  >  F  f  4 .  2  X .  4  H  2 1  =  , 

1  E10.4i2Xit.HRN  SO  =  ,  E  1  0  .  4 , 2  X  ,  7  H  A  L  P  H  A=  ,  F  4  2  ,  2  X  ,  «■  H  B  E  T  A  =  ,F4  2  .- /  /  > 

READ  <  1 0 '  1  )  P 
C 

C  INITIALIZATION  OF  TEMPERATURES 
C 

DO  5  I  =  1  , MM  AX 
DO  5  J  =  i  . N  M  ft  X 
W  1  (  I  .  J  )  =  A  L  ?  H  A 
U  2  <  I  i  J  )  =  3  E  T  A 
5  CONTINUE 
/  T  K  R  =  T  K  S  /  T  K  G 
V  M  I  =  ft  M  ft  X  -  1 
N  1  =  N  li  ft  X  -  1 
D  X  =  1  .  /  F  L  0  ft  T  <  M  1  ) 

D  R  =  1  /  F  L  0  A  T  <  N  I  J 

J  0  ft  =  0 

G  =  F  L  0  ft  T  (  ,M  I  -  1  >  t  FL  0  A  T  (  N  !  -  1  ) 

C 

C  DETERMINATION  OF  COEFFICIENTS  IN  DISCRETIZED  EQUATIONS 
C 

7  J99=J99+1 
DT I =0 . 

DT  2  =  0 . 

D  E  R  1  =  0  . 

D  E  R2  =0 . 

E  X  A  =  0  . 

E  R  A=  0 
C 

DO  10  I =2  .  M  1 
I  P  1  =  I  *  1 
I  M  1  =  I  -  1 
1  P2= 1  +  2 
I M2=  I  -2 
DO  2  0  J  =2  .  N  I 
J  P  I  =  J  +  I 
J  li  I  =  0  -  1 
J  P  2  =  J  +  2 
J  N  2  =  0  -  2 
C- 


-J*. 


c 

c. 

c. 


CALL  S  Q  U  E  L 
R  =  F  L  0  A  T  •:  0  -  1  )  <  0  R 


0  S  *  <  P  <  I  ,  J  P  1  >  -  ?  <  I  i  0  M  1  ?  )  /  \  2  .  *  D  R  *  R  t 

V  S  =  -  (  P  <  I P  1  ,  J  )  -  P  C  I  M  1  .  J  )  )  *  R  T  /  (  2  *  D  X  *  R  *  H  T  ) 

V  R  =  $  Q  R  T  <  LISHlSt-  VS  *  VS  ) 

R  £  B  =  R  t  *  D  ?  *  V  R  /  P  T 

C 

C  THE  DETERM!  N  =  T  I  OH  0?  THE  EDDY  D  I  F  F  0  $  I  V  I  T  Y  ,  E  ,  AND  1 
C  USING  A  MIXING  LENGTH  MODEL 
C 


V  $  U  =  0  . 

V  $1  =  0. 

0  $  L  =  U  S 
U  $  R  =  0  . 

IP  <  I  .Ed.  Ml)  GO  TO  50 

V  $  U  =  -  R  T  *  <  P  <  I  F  2  »  J  >-?<  1  »  J  )  )/<  2.  *DX*R*HT  ) 

5  0  IF  ■:  I  .  E  3  .  2  )  GO  TO  5  2 

V  $  L  =  -  P  T  *  ■:  P  :  I  ,  J  )  -  P  (  1  M  2  .  J  )  >  /  <  2  .  *  0  X  *  P  *  H  T  1 

5  2  IF  -I  J  .  E  3  .  2  >  GO  TO  5  4 

USL=(P(  I  .  0  )  -  P<  I  ,  0  M  2  .)  )  /  (  2  .  +  D  R  *  R  ) 

5  4  IF  <  J  .  E  5  Ml)  GO  TO  5  6 

(0  S  R  =  <  P  <  I  i  J  °  2  )  -  P  v  I  >  J  )  )  /  <  2  .  *  D  R  *  R.  > 

i  4  CONTINUE 

D  V  s  X  «<  VS  -J  -  V  3  L  >  /  •:  2  .  ■*  D  X  ) 

DU?R  =  :  USR-OSL  )/<  2  .  +  I) P  ) 

D  D  V  S  X  =  < V3U-2  .  » VS+VSL )/( DX*DX  > 

D  D  U  ?  R  =  <  U3R-2  .  *  US+USL  )/<  DR*DR  > 

22  =  1. 

23  =  1  . 

IF  <  D  V  $  X  .LT.  0  )  2  2  =  -  1  . 

IF  <  D  J  S  R  . L  T  .  0.)  Z  3  =  -  1  . 

C 

5  X  -Z  1  *0P  f  DP*  PS  *  V  J  s*  ftSS  <  D  VOX  >/  <  R  T  *HT  ) 
c  R  =  2  1  *  D  P  *  D  D  *  p  E  *  V  I  S  *  A  S  S  C  D  U  S  P  5  /  <  R  T  *  P  T  ) 

E  X  X=  2  1  *  ?.  2  *D  P  *  0  P*  P  E  *  V  I  3  *0  D  V  S  X/  <  9.  T  *H  T  > 

ERR  =  21  *23  *DP* OP* RE *V! 3  *  0  D  0  S  R  /  <  R  T  *RT  ) 

E  X  ft  =  E  X  A  +  E  X  7  G 
ERft=ERA+ER'C 
C 

C  FIND  N  U  ■:  R  E  .  E  :  R  )  > 

C 

C  =  1  .  -POR 
0  2  r  0  *  *  0  333  2  3 

01=2. -3. +02+3  *02**5  -2  *02**4 

PNI>=1  .<■*<  (REP  ♦  PR/  (POP*  0!  >)*<  1  -0**1  .  00467  )>+*0 
C 

C  CALCULATE  COEFFICIENTS  IN  EQUATIONS 
C  ■ 


A  A  1  =  RT  *<  POR  +  PP*EXXV  I  S  ,V<  RE  *PP*HT  *HT  ) 

A  A  2  =  R  T  *  <  1  .  -  POR  >  *  T  K  R  /  <  RE*PR*HT"HT  ) 

CC 1 =RT  *( POS*PR*ER/V IS  )/<  RE  *PP*RT  *RT  > 

CC2  =  RT  *  <  1  .  -  POR  )  *  T  K  R  /  C  RE*PP*RT  *  R  T  .) 

H  =  6  *  <  1  .  -POR  )*PNO*RT /  \ RE*PR*Dr+DP  ) 

3  =  RT  *  0  P  *  <  1  .  -  P  M  A  X  A  0  )  /  <  RE*PR*TF.  G*T  0  ) 

VI  =  -CCl  /  R-RT  P  0  R  R  +  P  R  *  E  R  R  /  V  IS  )/(  RE  *  P  R  *R  T  *  P.  T  )  +  V  S 
V2=-CC2/R+RT  *PORR*TKR/< RE*PR*RT*RT  ) 
,J1=US/HT-RT*(P0R>:*PR*EXX/VIS>/<RE*PR*HT*HT) 
U2=RT  *TKR*PORX/\  R  E  *  PR  *  H  T  *HT  ) 

C 

C  COEFFICIENTS  USED  IN  UPWIND  DIFFERENCING 
C 

IF  <  U 1  . LT .  0 .  )  GO  TO  00 

A  8  =  -  1  . 

3  8  =  1. 

C  8  =0  . 

GO  TO  43 
6  0  A  8  =  0  . 

8  ft  =  -  1  . 

C  3  =  1  . 

63  IF  <  VI  . LT .  0.  >  GO  TO  65 

A  8  =  -  I  . 

88=1. 

C  9  =  0 

GO  TO  67 
65  A  9  =0  . 

8  9  =  -  1  . 


D  E  R  1  V  A  T  I  V 


33333 


/RT 


■  I  l  - 


C  A  =  ! 

&  ?  CONTINUE 

c 

A  t  t  J  >  -  -  C  C  1  /  C'  R  -  D  P  +  V  1  *  A  ?  /  D  R 

9  1  (  J  )  -  2  *<  4  A  J  /  D  X  .• ■'  D  X  *  C  C  1  /  t>  R  /  D  ft  >  +  H  «•  0  1  *  B  3  /  D  X  ♦  V  I  *  B  ?  /  D  R 

C  I  (  J  )  -  -  c.  C  1  /  D  R  D  P  +  V  I  *  C  9  /  0  R 
DH  J  )  -  H 

E  i  (  j  >  =  A  a  i  •  y  :  ■  i  ?  l  ,  j  )  +  y  i  <  i.n  l  .  J  >  >/  ox  /ox  -u  i  ft  3  *y  i  <  i  m  l  .  J  )  +cs  *w  :  ■:  I  p  l 
1  J  )  >  /  D  X 

c 

C  USE  CENTRAL  DIFFERENCING  FOR  THETA  X  AMO  THETA  R 

C 

98=0 
9  a  =  n 
C  8  =  6  5 
C9  =  0  .  5 
A  8  =  -  0  .  5 
A9  =  -0  .  5 
r 

A  2  (  J  >  -  -  C  C  2  /  D  P  /  V  R  +  V  2  *  A  9  /  D  R 

3  2  (  j  )  =  2  *"  <  AA2/DX/DX+CC2/DR/DR  >  +  H+  U2*B8/DX  +  V2*83/DR 

r  2  <  J  T  =  -  C  C  2 /  D  R  DR  ♦  V2  *C9  / 0  R 
0  2  (  J  j  =  H 

E  2  (  J  )  =  A  fi  2  *  ■;  y  2  (  I  P  1  ,  j  >  +  U  2  C  I  K  1  ,0  >  )  /  D  X  /  0  ",  -  U  2  *  A  3  *  y  2  (  IK!  ,  J  1  +  2  3  *  W  2  •:  IP! 
1  J))/DX+8 

CO  CONTI H U E 

C  ADJUST  COEFFICIENTS  SO  AS  TO  SATISFY  BOUNDARY  CONDITIONS 

B  I  <  2  5  =  8  1  •!  2  )  ♦  A  1  (  2  > 

S  2  C  2  D  =  e  2  •:  2  >  +  A2(  2  > 

B  1  <  N  1  )  *  B  I  <  N  1  )  +  C  1  (  N  1  )/  <  1  .  ♦  D  R  *  R  T  +  R  N  G  W  ) 

E  1  (  N  1  )  =  E  !  (  S  1  )  -  C !  <  N  !  )  *  D  R J'  R  T  *  R  N  G  U  *  7  A  /  <  T  0  *  !  !  .  «-  0  P  *  R  T  *  R  M  G  W  J  > 

8  2 ( N 1  ) =B2<  SI  >+  C2 ( N 1  )/<  1  .  ♦OR*RT*RMSW  ) 

S  2 ( N i  >  *  E  2  t  N 1  )  -  C  2  < HI  >  D  R  *  R  T  *  R  N  SO  *  T  A  /  <  T  0* t  1  .  +DR*RT*RNSW)  ) 

C 

C  EXECUTE  SOLUTION  TECHNIQUE  TO  SOLVE  COUPLED  TRIDI AGONAL  SYSTEM 
C 

CALL  TRISBU  (  N  (1 A  X  .  S  OL  N  1  ,  SO  L  M2  > 

C 

C  EXECUTE  RELAXATION  TECHNIQUE  IN  ORDER  TO  UPDATE  UNKNOWN  VALUES 
C 

DO  3  0  J  =  2  »  N 1 
C. 

s  n  l  =  y  i  <  i ,  j  >  *  us  *<  s  oi  n  i  <  o  >  -a  i  <  i  >  j  )  :> 

S N  2  =  W  2  (  I , J  >  +  U$*< SDL  N2< J  )-V2C I  > J  )  > 

C 

ERR1  =  A  8  S  -I  S  N  I  -  S  0  L  N  1  <  J  >  0 
E  F:  R  2  =  A  8  S  <  S  N  2  -  S  0  L  N  2  <  J  )  > 

D  T  1  =  D  T  I  +  E  R  ?.  1  /  G 
/  DT  2  =  DT  2+ERR2/G 
C  V 

C  FIND  MAXIMUM  ERROR  IN  THE  NET  AND  ITS' POSITION 
C 

CALL  C0MV2-:  I  ,  J  ,  £  R  R  1  ,  E  R  R  2  ,  L  I  1  .  L  1  2  >  L  J  1  .L02.DERI.DER2) 

C 

9  1  (  1  »  J  0  =  S  N  1 
y  2  <  I  ,  J  )  =  S  N  2 
30  CONTINUE 

v  i  ( i  .  i  o  =  a  i  •:  i  >  2  > 

a  2  <  i  >  l  >  =  a  2  •;  i  j  2  > 

ai(I  >  N  N  A  X )  =  (  a  1 < I  , N1  )  +  DR*RT  *RHGW*TA/T  0  )  / <  1  .  R * R T  + R N G » > 

U2(  I  i.NMAX  >  =  <  W2(  I  <  N  1  J  +  DR’»RT*RWSW*TA/TO>/<  1  .  +DR*RT*RNSU) 

10  CONTINUE 

DO  3  5  J  =  1  >  N  M  A  X 
y  1  <  M  M  A  X  >  J  )  =  W  I  (  M  1  ,  J  .> 
a  2  (  M  M  A  X  ,  J  0  =  a  2 1  M  1  .  0  > 

a  i  <  i ,  j  )  =  i  . 

V  2(  1  ,  J  >  =  •:  DX  *HT  *8N  SO  +  y  2  (  2  »  J  >  >/(  J  .  +D  X  ♦  H  T  +  R  N  S  0  0 
35  CONTINUE 

IF  •:  J  3  8  .  L  E  3  0)  GO  TO  33 

GO  TO  39 

33  J  RITE  (  103.300)  J99.LM.LJ1.DER1.LI2.LJ2.DER2 
300  FORMAT!  1 X  »  2  H  K - »  13.  3H  X= .  12,  3H  P= .  12.  ?H  ERRHBXI*.  E12.6, 

!  2  .X  .  2  A X  =  .  I  2  ,  3  H  R  =  ,  12.  ?H  ERFM  AX2  =  ,  E  I  2  .  S  .  /  ) 

a  R  I  T  E  :  I  0  3 . 4  0  0  >  D  T  i  ,  D  T  2 

400  FORMAT-:  IX  .  3HERRAVG1  =  ,  E  I  2  .  e  -  !  OH  ERRAVG2-,  EI2  &./) 

3Q  CONTINUE 

IF  <vT?9  .GT.  J  33  M  AX  )  GO  TO  40 


IF  (Ml  .  G  T  .  EPS)  GO  TO  7 
IF  •:  P  T  2  .  G  T  EPS)  GO  TO  7 
4  0  C  OS T I N U  E 

IF  <.!??  .  LE  .  ?0>  GO  TO  45 

WRITE!  10  3.300.'  J  ?  ?  j  LI  1  ,  L  J  1  .CEP  1  ,LI2,LJ2»D£R2 
■J  R  I  T  E  I  1  0  3 . 4  0  C  )  0  T  1  ,  D  T  2 
45  CONTINUE 

£F  =  33  RT  £SA*E>:fl*ERft*ERfl> 

WRITE  (  108,  425  )  F  F 

425  FOR OAT  ( IX. 13HEFFECTIYE  E -  ,  E  1  0  .  4  ,  /  ) 

C 

C  THE  SOLUTION  FOLLOW? 

C 

WRITE  (108,450) 

450  FORMAT* /  , IX, 47HTEMPERATURE  DISTRIBUTIONS  FOR  THE  GAS  AND 
DO  55  K=l»MMftX 
L  =  MM  fiX  -  K  +  1 

WRITE  (108,500)  (  U  1  (  L  .  J  )  >  0  =  I  ,  NM  A  X  ) 

WRITE  (108,300)  (W2(L,J)»J=l,NHflX> 

500  FORMAT  (  2X,  1  1  r  9  .  4  ,  /  ) 

800  FORMAT  (IX,  1  H  +  ,  1  1  F  9  .  4  ,  /  /  ) 

55  CONTINUE 
C 

WRITE  (10*3)  LI  1 
WRITE  (10*4)  « 2 


*  * 

FR, CL,  SL . UL  ,( HOTBSV  ) 


,5  03 


IPLANTATI  ON 


0000 
i  A3  4 
38 
bOAC 
5  DC  A 


» 

c 


5  E  2  2 
5  E  5  6 
5  EGO 
5  E  7  2 
5  E  3  2 
5  E  C  2 
5  F  3  4 
5  F  3  S 


5  F  3  2 
5  F  A  4 
659  8 
659C 


6  5  A  2 
6  53  H 
864  4 
6  64  8 


6  65  4 
6658 
6  65  E 
6  6  3  2 
6  6  3  6 
6  6  3C 
6  69  E 
6  72  4 
©  8  <;•  A 
6  8  3  E 


Dc  S  SEGMENTS*** 


SOLID,//) 


DD  = 

0  .  0  0  50000 

RT  =  0  . 

0300000 

HT  * 

0  1000000 

DSNS 

- 

1 2 1 1 E*0 1 

v i sco= 

.  1  5  3  0  E 

-04 

RET 

*  40  00 

GX  = 

0 

02  50  0 

SK  = 

0 .  7  5  OP  C 

0*  3 

3  ?? 

TEX  = 

20 

o 

» — 

o 

*15.00 

H  = 

31 

H  = 

3  1 

NF  = 

1  CO 

N  G  = 

2V 

OVER 

RE 

L  AX  AT  ! ON 

F ACT  OR  = 

1.05 

TEST* 

0 

0  1 

NO  U  G 

= 

Q  VU  * 

0  .  00 

MU  41 S 

= 

«***«*» 

q  y  s  = 

0  .  00 

B  l  0  X  A  G  E  1  ♦ 

P  H  I  N  A  0  »  010  0  P  9  A  X  A  0  *  0  100 

11100  =  11  12200  =  21  J11A0  =  11  J4A0  =  21 

B  L  0  X  A  G  E  2* 

PfllNAl  =  0  1 00  PHAXA1  =  0.100 

11101  =  50  12201  =  50  01101  =  50  J4AI  =  50 


PR* 


0.72 


J«ON«X--  200  EPS*  0  000  1  2  1  =  1427E*D0  RNSO=  OOOOEOO  ALPHA*  1.00  BE  T  0  = 


K  =  1  X  =  3  0  R  =  2  6  E  R  R  M  A  X  l  =  . 2  6 1 2 1  1 E -  0 1  X  =  30 
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